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Abstract: 


The mathematical modelling of electric current flow in excitable tissues is briefly reviewed and three 
fundamental problems of the stimulation by external electric fields are defined and posed, in terms of 
suitably space-averaged fields. Then, it is given an outline of a new analytical approach to study just- 
threshold dynamics when certain target elements in excitable tissues are stimulated by external electrodes 
or magnetic coils. Both time and spatial aspects of membrane’s non uniform polarization are considered. 


In principle, the proposed approach can be applied using any electric membrane model (not restricted to a 
capacitance in parallel with a non-linear conductance). The non-linear ionic current of the corresponding 
membrane model must be approximated by a suitable multivariable polynomial (expressed in the state 
variables of the membrane: channel activation variables, membrane voltage, channel recuperation 
variables and accommodation variables). 


The concept of “region of influence” of an external field over the target element is introduced as a key 
concept that allows a nonlinear modal analysis of membrane’s polarization in the target element. In the 
fiber’s case we have a one-dimensional interval of influence. In the electric syncytium case, we have a 
two- or three-dimensional region of influence. 


Membrane voltage, activation variables and recovery variables are assumed to remain in their rest values 
on the boundary of the region of influence until a threshold pattern is achieved, inside the region, due to 
the perturbation produced by the applied electric field. 

So, the nonlinear partial differential equations that describe the stimulation of the target element (fiber or 
electric syncytium) can be studied inside the region of influence, with homogeneous boundary conditions, 
relative to the constant rest state of the membranes. 


Then modal analysis allows the representation of the different fields (functions of space and time) like 
membrane voltage, activation variables, recovery variables, and so on, as linear combinations of known 
space functions with unknown time dependent mode amplitudes: the modal representation of the fields 
inside the region of influence. The known space functions belong to the complete set of orthogonal 
eigenfunctions of a suitable posed Sturm-Liouville problem. 

Substituting the modal representation of the fields in the nonlinear partial differential equations that 
describe the stimulation process, a set of nonlinear ordinary differential equations in the mode amplitudes 
is derived. Each one of these ordinary differential equations has an additive term that depends of the 
effect of the external electric field (produced by the active electrode or magnetic coil) on the 
corresponding mode amplitude. 


The excitation dynamics due to a perturbation caused by an imposed external electric field can be studied 
as a trajectory in the space of mode amplitudes, from the rest state and up to the threshold states beyond 
which an action potential emerges. In a suitable truncated mode amplitude space, a threshold manifold, a 
decaying set, and an amplifying set are introduced to study the emergence of an action potential. 


By means of this theoretical framework it is possible to pose excitation problems as control problems in a 
suitable state space and the introduction of optimal control criteria (such as least energy dissipation in 
tissues or in the whole stimulating equipment-tissues system). 


The threshold value of a perturbation (from the rest state) that causes instability in the space of mode 
amplitudes is estimated using a modified version of Eckhaus method. 
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The classical one and two factors excitation theories are derived, and the phenomenological “factors” of 
these theories are reduced to modal amplitudes. 


It is then outlined a general procedure to derive strength-duration formulae, giving the derivation of 
Lapicque-Hill’s formula, Lapicque-Weiss’ formula and a new formula that could be used in the current 
practice of cardiac pacing. 


An extension of the Bonhoeffer-van der Pol-Fitzhugh model for an excitable membrane is proposed, and 
then applied to the modal analysis of the resulting nonlinear cable equations for excitable fibers and 
excitable electric syncytia. A “three factor theory” of excitation (electric or magnetic) is derived. 


Also, it is pointed how a solution could be given to several old problems of electro stimulation (related 
directly or indirectly with strength-duration curves) that have remained open until now, using the 
theoretical framework proposed in this paper. 


The concept called “family of threshold states” (for each given target element) is proposed and its 
applications to synaptic excitation and to spontaneous electrical activity of excitable tissues are suggested. 


The main shortcoming of the proposed general procedure is that it is unable to cope with the growth and 
propagation of the action potential, being restricted to just-threshold dynamics. 

However, from one side this may be enough for the discussion of differential excitation and related 
subjects of interest from the viewpoint of clinical electrophysiology, while from the other side the 
obtained results could be used as guidelines for doing new experiments and new studies of digital 
simulation of threshold dynamics. 
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1. Introduction 


Strength duration curves measure the excitability of a certain target element in each 
electrode-tissue or stimulating coil-tissue system. 

The target element may be a single cell (such as a nerve or skeletal muscle fibre) or 
even a whole network of electrically coupled cells (such as the myocardium or visceral 
smooth muscle). 

The response of the target element is all or none emergence of a propagated full action 
potential. 


When determining strength-duration curves in electric stimulation, the stimulus is a 
rectangular pulse of electric current of amplitude I and duration d, applied to the volume 
conductor formed by the biological tissues, through an array of external electrodes. 

The resulting field of electric current density polarizes non-uniformly the excitable 
membranes of the target element, so that the classical concept of a uniform membrane 
threshold potential cannot be applied to study the emergence of an action potential. 

In synaptic excitation and in excitation by intracellular microelectrodes -that are both 
essentially localized and spatially non-uniform phenomenon- a small current source is 
localized at the membrane level or in the intracellular space respectively and produces 
either a depolarization or a hyper polarization of the whole membrane. 

But in the excitation by external electrodes there aren’t sources or sinks of electric 
current in the bulk of the tissues, where the target element is located. 

Because of this, part of the membranes of the target cell (or cell network) will be 
depolarized and the other part will be hyperpolarized. 


If the target element is a fibre stimulated by an external active electrode that is close 
enough to the fibre, besides the well-known cathode make and anode break excitation, 
two new ways of excitation appears, related to the always dual polarization of fibre’s 
membranes: cathode break and anode make excitation (1). 

Cathode break excitation presupposes a blockade of the previous cathodic make 
excitation, a phenomenon known as “anodic block from the cathode” (2). 

Each has its own strength-duration curve. 


From an experimental point of view the strength-duration curves may be obtained fixing 
the duration and varying the amplitude of the applied current pulse until threshold 
amplitude I, is found. Doing the same thing for different pulse durations, a relation 
between I, and d can be established. 

When I = L,, the excitation will succeed in half of the cases. Varying I for a fixed d 
the percentage of successful excitations will pass from 0 to 100% in a small 
neighbourhood ofI,. So, under normal physiological conditions we can idealize what is 
a threshold fringe by a threshold curve (3). 

This curve may be summarize by a few phenomenological parameters: the rheobaseR, 
the limit threshold chargeQ,, the o parameter and the system’s time constant t, (see 
Fig.1) 


R= lim I, (dd) , Q.,=limd.I, (dd) , t 
t t alo t 


Sok, 
d>+0 sR R 
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Fig.1 Strength-duration curve in log-log form 


Instead oft, an older time measure is still widely used in functional electric stimulation 
of nerve and muscle, and in cardiac pacing: the chronaxie, defined as the pulse duration 


( 
t., such that (1, 3): Sag ne 


In the same electrode-tissue system we may have several target elements, so that we can 
construct a strength-duration curve for each one, and these curves will be, in general, 
different from each other. 


But if we change the position or type of electrode, leaving the rest of the things as they 
were before, the new strength-duration curve for a given target element will be, in 
general, different from the old one. 

Indeed, the results of the experiments show that the parameters of the strength-duration 
curve depend not only on the excitability properties of the membranes of the target 
element but depend also on certain geometric and bioelectric parameters of the system 
formed by the electrodes, the tissues, and the target element. For example: radius of the 
excitable fibre, distance between the electrode and the target element, size, and shape of 
the active electrode (if it is located close enough to the target element), amongst others. 
Even the chronaxie, introduced by Lapicque as an intrinsic time measure of the 
excitability of the tissues, depends also on the radius of the electrode, the distance 
between the electrode and the excitable tissues, and other parameters of the electrode- 
tissue system (4, 5, 6). 


So, strength duration curves measure something that may be called a “global 
excitability” of the target element in the framework of a given electrode-tissue system, 
by opposition to the “local excitability” of the membranes of the target element (7). 
This local excitability could be measured, at least in principle, determining the strength- 
duration curves for the uniformly polarized membranes (by a patch-clamp technique, for 
example) and relating its parameters with the known properties of the ionic channels. 

A full explanation of the experimental results that show a relation between the 
phenomenological parameters of the strength-duration curve and certain parameters of 
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the electrode-tissue system has remained as an open problem from more than seventy 
years. 


This problem has more than an academic interest, because it is closely related with two 
basic questions of clinical electrostimulation (8, 9): 


-Which elements will be excited by a given current pulse in each electrode-tissue 
system? 

-How to reach a selective stimulation of certain target elements in a whole excitable 
tissue? 


To get a solution to these two problems, it seems necessary to solve first the following 
three problems of the electric stimulation by external electrodes: 


(1) Given a certain electrode array in a certain volume conductor formed by 
biological tissues, to find the resulting field of electric current density. 


(2) To obtain a relation between this current density field and the space-time 
pattern of polarization of excitable membranes, including the field of 
transmembrane voltage Vu and the whole set of activation, recovery, and 


adaptation variables of the ionic channels, represented by the vector field W 


(3) To identify certain spatial patterns of ( Vu,W) that in some sense could be 
considered as “threshold states” of the excitable membranes in the target 
elements. 


If we could solve these problems using an analytical approach, then we could expect to 
derive some meaningful formulae for the strength-duration curves. 

The aim of this work is to give a brief outline of a general analytical procedure to solve 
the three problems of electrostimulation and to do the derivation, for a single fiber and 
for a whole network of electrically coupled cells. 

Two mostly empirical formulae are used to characterize the global excitability of nerve 
and muscle during functional electric stimulation, and the excitability of myocardium in 
cardiac pacing: 


I 1 
The classic Lapicque -Hill’s formula (7) = = a 
es 
. . . t t. 
The also classic Lapicque-Weiss formula (6) ie 1 + a 


We derive and discuss both formulae to enlighten the general procedure to be followed. 
As a further consequence of this derivation, a new formula for the strength-duration 
curve during cathodic make excitation of ventricular myocardium is obtained. 
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2. The electric current flow produced by external electrodes and magnetic coils in 
excitable tissues. 


Let us first consider the electric current flow produced by the external electrodes. If the 
emergence of an action potential were a short-range spatial phenomenon, we would 
need a detailed description of the flow in the environment and in the interior of the 
excitable cells. But fortunately, this emergence is a long-range spatial phenomenon 
under common physiological conditions. So, it is possible to work with suitable space 
averaged fields, smoothing unnecessary details and greatly simplifying the 
mathematical formulation of the three problems of electrostimulation using a 
reformulation of the so-called bidomain model. 

The bidomain model is a phenomenological approximation (10, 11) originally devised 
to describe the electric current flow in cardiac tissue in the framework of 
electrocardiography and magnetocardiography, but it can be restated to study threshold 
dynamics in any electrode-tissue system under (12, 13). It assumes that the discrete 
nature of the tissues can be ignored in favour of a simplified continuum approach. Three 
continua are postulated in the original version of the theory, each one filling the same 
spatial region: an intracellular domain, an interstitial domain, and the continuum of 
excitable membranes. In each continuum, several fields are introduced that verify a 
coupled set of nonlinear partial differential equations. All this can be suitable founded in 
the equations of electrodynamics and electrolytic transport, taking into due account the 
restrictions stemming from the complex geometry of the tissues, as has been recently 
shown (12). 

To study the problems of excitability it is convenient to substitute the intracellular and 
the interstitial domains by a single ohmic volume conductor, the equivalent volume 
conductor. As this new continuum is independent from the continuum of excitable 
membranes, it is possible to pose and solve the problem of finding the electric current 
distribution produced by an array of external electrodes as in any other ohmic volume 
conductor. 


Under quasi steady conditions the field of electric current density J verifies the 
equation 


V.J=0 [1] 


(charge conservation). If ® is the field of electric potential and G is a suitable tensor 
field of local conductivities, 


J=-G.VO [2] 
From [1] and [2] follows: 
V.(G.V®)=0 [3] 


This equation may be applied in the region B occupied by the volume conductor. 
(Here. represents the scalar product and V.J is the divergence of the vector field J and 
V ® is the gradient of the scalar field ®). 


Adding boundary conditions (for example, assuming that the component of J normal 
to the boundary of the region B vanishes at points not in the surface of an electrode and 
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that either this component or the value of ® are given at each point in the surface of 
each electrode) the scalar field ® can be found by well known analytical or numerical 
methods. 
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3. The response of the target element to electric and magnetic stimulation 


Now we must turn to consider the target element. If it is a single excitable fiber, it can 
be represented by a global, spatially one-dimensional, and nonlinear cable model (14, 
15). The cable may be idealized by a certain curve in the region occupied by the volume 
conductor. If Ss represents the arc length along that curve, a general version of the 
nonlinear cable equations for a uniform fiber is the following (13): 


O Nu — Ooo Ne foo 
Pao + Ral WIR 4 pu Ds? [ 4a] 
OM =F Vy,W) [4b 
Here ty, = Ry Cu is membrane’s linear time constant, R,, is the resistance of the 


unit area of membrane, C,, is the capacitance per unit area, J:V,,.W) is the ionic current 
density through membrane’s ionic channels, 2 ,, is membrane’s space constant, 
F (V,,,W) represents the kinetics of the activation and recovery variables of the ionic 


channels and its scalar components are given by expressions of the Hodgkin and 
Huxley’s type (but that vary from one membrane to another), and t represents time. 
The coupling between the state variables ( V,,,W) of the nonlinear cable equations [4a], 


[4b] and the applied electric current field in the equivalent volume conductor is given by 
2: 


2 
the source term/ y” — where e is equivalent to the activating function 


ae Os? 
introduced by Rattay (2). Under quasi steady conditions 
2 
“2 =F) [ac] 


I(t) is the total applied current and F(s) is what we call the form factor (13, 16). When 
the fiber is bounded, another coupling with the external current field appears through 
the boundary conditions at the end points of the fibre (13). 


Instead of a point electrode close enough to the target fibre now we consider the case of 
the so called “fluid electrodes” that produces a uniform field of electric current density 


in the region in which the target element is located. Then the field - V M(t, r) = E,(t) 
is space independent. 


> > 
If r = r,(s) is the vector function that gives us the position of the (perhaps curved) 
fiber in the equivalent ohmic conductor (which is assumed to be isotropic and 


d > > > d- 
homogeneous), we obtain ao = - E,(t). t,(s) , being t,(s) = ast the 


ae 
targent unit vector to the fiber considered as a space curve. If n,(S) is the unit normal 
vector and K,(s) is the curvature of the fiber in the considered point, the activating 
function is given by 

d? 


ds? 


(s) =-K, (s) Ep(t)-n,(s) [4d] 
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(We have employed the first of Frenet’s formulae (21) to derive the preceding 
equation). So, if the external field is not parallel to the fibre at the considered point and 
if the fibre has a non-zero curvature there, the activating function will be different from 
zero. But even if the fibre is straight (so that K,(s) = 0 everywhere and the activating 
function is always zero), the fibre will be polarized by the applied field due to the 
boundary conditions at the ends of the fibre. Of course, this presupposes that the fibre is 
finite. Please note that the stimulation of target elements located deep in the central 
nervous system could be studied using this framework and the region of influence 
concept, introduced later, in part D. Also note that some still apparently unexplained 
results obtained by Rushton in 1927 (17) with curved fibers stimulated by fluid 
electrodes could be tackled analytically for the first time using the nonlinear modal 
approach proposed in this paper. 


Of course, digital simulations applying suitable computer packages can always be done 
(and as a matter of fact has been done and are done), to cope with the same problems 
from a numerical point of view, for some selected values of the system parameters. 

The powerful software for plotting multidimensional data and examining the various 
possible projections might enable the user to extract or suspect algebraic results form 
numerical information, although this is not always feasible. 

Nevertheless, analytical models, if possible, allow a deeper understanding of the 
interrelations in the modelled system and are very useful in the design of better digital 
simulation runs. 


If the target element is an electric syncytium, it can be represented by the equations for 
the continuum of excitable membranes. If, as in ventricular myocardium, the anisotropy 
is unequal (i.e., the anisotropy ratios of the interstitial and intracellular spaces are 
different) the nonlinear membrane equations are coupled to the electric current field in 
the equivalent volume conductor by a source term that is, like Rattay’s activating 
function proportional to I(t), but much more complex. This term vanishes if the 
anisotropy is equal and uniform. As in the single fiber case, a coupling between the 
equivalent volume conductor and the excitable membranes always appears through the 
boundary conditions, now at the surface of the continuum of excitable membranes. The 
boundary conditions are of the Neumann’s type, proportional to the total applied 
current, so they are nonhomogeneous while the external current is flowing through the 
equivalent volume conductor, and homogeneous when this current vanishes. 

The equations for the continuum of excitable membranes are the following (12): 


Ln (Ty +R, ICV, W)) = VR CV W4¥| MEG, 1.6090] [5a] 
OW =. = 
“_™ =F(V,,,W 5b 
A (Vu. W) [5b | 


Here Ym is the area of excitable membranes per unit volume, f, and f; are the volume 


fractions of interstitial and intracellular domains respectively, G, and G, are the 


e 


Ry 


(tensor) conductivities of these domains, is a tensor space constant and the rest 


m 


of the symbols have the same meaning as in the single fibre’s case. 


10 
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Moreover, C =f,G,.G'.f,G, where G~'is the inverse of the operator 
C28 E086 10) 13), 


The boundary conditions are the following (assuming that in bidomain’ s boundary the 
electric current goes to or comes from the surrounding tissues only through the 
interstitial domain) (12, 13): 


~ 2 ~ ad 


n.C.VV,,=n.(, G,.G7.J) [5c] 


= 
where J is the electric current density in the equivalent ohmic conductor, produced by 


x 
the external electrodes, and n is the outer normal vector at the corresponding point in 
bidomain’ s boundary. 


Before applying these equations to the polarization of excitable membranes, it might be 
appropriate to briefly review how they are derived from standard bidomain equations. 

In bidomain’ s model it is assumed that the electric current densities in the intracellular 
and interstitial domains are given by the generalized Ohm’s laws: 


J, =-f,G,.VV. [sd] 


J, =-f,G,.VV. [Se] 
In these relationships, V, and V, are suitable electric potentials. Equations [5d] and [Se] 
are coupled through the relations. 


Vii 7d, [sf] 


VJ. =- 4m Ju [5g] 


These relations involve the total electric current density through the continuum 
excitable membranes: 


O 


a Vu tI W) [sh] 


Ju =Cur 


As the operators G, and G, are symmetric and positive definite, they are invertible, and 
they can be diagonalized. Assuming that the principal directions are the same both for 
extracellular and for interstitial domains, these and other operators constructed from 
these ones by algebraic operations) will commute. 

Up to this point we have followed the standard approach (10). Now, we retain the 
continuum of excitable membranes, but we substitute the interstitial and intracellular 
domains by a new continuum: the equivalent ohmic volume conductor, as follows. We 
introduce a new potential @® such that it verifies the equation (being 
G = £,G, + f, G,): 


G.VO=f,G,.VV, +f, G,.VV [5i] 


11 
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(For details, see (12)). By construction, the field ® verifies equation [3], because 


3 
J; + J, is the total electric current density in the tissues. Then, considering that 
Vu =V;-V [si] 


is the field of membrane voltage, we can obtain V V, and V V, as functions of V V,, 
and V ® from equations [5i] and [5j]. Then, substituting the thus obtained formulae in 
equations [Sf] and [5g], subtracting equation [5g] from equation [5f], we derive, after 
some transformations (including a multiplication of both members by R,, ), the equation 


[5a] as desired. 


Now we can turn to see how equations [Sa], [5b] and [5c] can be applied to the 
polarization and excitation of a functional syncytium, such as cardiac muscle or 
electrically coupled smooth muscle. 

Let us consider, for example, cathodic make excitation, begining with the membranes of 


the target element in their rest state given by the constant values V, and We of the state 


variables. In principle we can solve the corresponding initial-boundary value problem, 
either for a single fiber (using equations [4] ) or for the continuum of excitable 


membranes (now using equations [5] ), obtaining the spatial distribution of V,, and W 
at the end of the current pulse, thus solving the second problem of electrostimulation. 
The rest state ( V,,We) turns to be a locally stable equilibrium solution of the nonlinear 


membrane equations when the external current vanishes. In principle also, we could 
establish the amplitudes of the perturbation that somewhat produce instability and allow 
the emergence of an action potential, identifying certain spatial distributions V,, and W 
that in some sense could be considered as just-threshold distributions, thus solving the 
third problem. Fixing certain duration d for a rectangular pulse and finding the 
amplitude I, that at the end of the pulse produces a just-threshold distribution of the 
state variables of the membranes of the target element would allow us to derive a 
theoretical strength-duration relation. 

Now, this direct approach doesn’t seem to be possible from an analytical point of view, 
even for the simplest case of a uniform excitable fibre. Linearizing the cable equation as 
was done by Rubinstein and Spelman (18), it is possible calculate approximate 
polarization patterns along the axis of the fibre. But with the linearization the threshold 
disappears! 

Some very interesting results have been obtained working numerically from the very 
beginning with the nonlinear equation, for the case of an unmyelinated nerve fibre (2). 
However, it seems difficult to continue the numerical ab-initio computations without the 
guidelines of a theoretical framework, such as is afforded by an analytical theory. 


12 
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4. Modal analysis of threshold dynamics: a brief outline 


If we want to continue with the intended analytical approach to threshold dynamics, we 
need a new idea that allows us to drastically simplify the original nonlinear 
mathematical problem. This idea happens to be the region of influence of the applied 
current field over the membranes of the target element. 

This region of influence can be suitably defined as a bounded and closed region such 
that outside of it, the polarizing effect of the external current field can be neglected. So, 
when a point of the target element is in the boundary of the region of influence, its 
membrane voltage and its ionic channel variables remain fixed at their rest values. 

If the target element is a single fibre and if there is a single active electrode (cathode or 
anode) close enough to the fibre, this region is simply the interval of influence of the 
electrode over the fibre (13, 16, 19). If the target element is a functional syncytium and 
if there is a single active electrode, like in cardiac pacing, we have a whole three- 
dimensional region of influence of the electrode over the tissue (13). 

Using the region of influence, it is possible to pose and then solve a spatial Sturm- 
Liouville problem (13, 16, 20) (with Dirichlet’s boundary conditions for a fiber and 
with mixed Dirichlet’s and Neumann’s boundary conditions for a bounded continuum 
of excitable membranes). 

The fields V,, and W can be developed as a series of functions in terms of the 
corresponding eigenfunctions or modes, with mode amplitudes that depend only of 
time. 

If we approximate membrane’s ionic current density J(V,,,W) and the scalar 


components of the function F(V,,,W)by suitable polynomials in terms of the 


deviations of the state variables from the rest state, the nonlinear partial differential 
equations [4] and [5] can be substituted by an infinite set of ordinary nonlinear 
differential equations in terms of the amplitudes of the modes. 

Moreover, it is possible to truncate these equations for the mode amplitudes from a 
certain mode index onwards. In that case the whole spatial distribution of the state 
variables V,, and W of the target element membranes can be represented by a single 
point in a suitable finite-dimensional space of mode amplitudes. 

The rest state is represented by the origin. Furthermore, the widely separated time scales 
of evolution of activation variables, of membrane potential and of recovery variables, 
allows the partition of that space in two sets: a decaying set and an amplifying set. The 
rest state is always in the decaying set. If at the end of an applied pulse the 
representative point is in the decaying set, the excitation will fail and the membranes of 
the target element will return to rest. But if it is in the amplifying set, an action potential 
will be produced. 

The amplifying set is separated from the decaying set by a more or less well defined 
transition region that can be idealized by a certain threshold manifold in the space of 
mode amplitudes (See Fig.2, Fig.17, Fig.24 and Fig.31 ). 


13 
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Amplifying Set 


Reachable Fs : 


Sets from 
the origin 


Pseudo-Threshold manifold 


Fig.2 Decaying Set, Amplifying Set, Pseudo-Threshold manifold and Reachable Sets in the space of 
modal amplitudes. 


Each point of this manifold corresponds to a threshold spatial distribution of the 
variables V,, and W, so that we obtain a definite answer to the third posed problem of 
the excitation by external electrodes. Moreover, to each history of applied current 
corresponds a certain orbit in the space of mode amplitudes (the effect of the external 
current on the rate equation for a given mode amplitude appears as a source term 
proportional to I(t)). This allows us to solve the second problem of electrostimulation. 
Considering the intersection between the decaying set and the region reachable from the 
origin by cathodic or anodic pulses, applied by a single active electrode, phenomena 
like cathodic make excitation, anodic make excitation, anodic block from the cathode, 
anodic break excitation and cathodic break excitation can be analysed by a thorough 
examination of the location of the reachable subset of the threshold manifold in mode 
amplitude space. 

A rectangular pulse of externally applied current will be just threshold if and only if, at 
the end of the pulse the point which represents the state of the membranes of the target 
element belongs to the threshold manifold. Now, linearizing the dynamics of the mode 
amplitudes in the decaying set, simplifying (if necessary) the threshold manifold, and 
considering the just-threshold condition for a rectangular pulse of applied electric 
current, it should be straight-forward to derive analytical expressions for the strength- 
duration curves, considering the effect of any number of activation and recovery 
variables as needed. 

Now, to be able to develop an analytical approach to the threshold manifold, it is 
necessary to reduce as much as possible the order of the dynamics. 

This reduction must be done in two steps. 

First, we must construct the region of influence with a proper size and shape. 

Second, we must try to identify some dominant mode amplitudes, such that the others 
amplitudes can be neglected as a first approximation. 

If the region of influence is too big, we will need too many eigenfunctions to reproduce 
the details of the spatial patterns of the membrane voltage and the ionic channel 
variables. 

But if it is too small, we will not be able to reproduce these spatial patterns, not matter 
what number of eigenfunctions we intend to use. 

A somewhat similar problem was faced by physicists at the beginning of nuclear reactor 
theory because neutron flux also vanishes gradually outside the core of a reactor. This 
problem was successfully solved introducing the idea of an “extrapolated distance”. 
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Using this concept, it is possible to define a virtual reactor core, greater than the actual 
core, and such that at the boundary of this virtual core the neutron flux can be 
considered as being always zero. 

Of course, in our case the underlying physics is quite different. Here we don’t have an 
equivalent to the core of the reactor. The spatial patterns of the fields don’t originate in 
the same kind of balance between internal production, absorption, and diffusion 
characteristic of the neutron case. But the always localized polarization of the excitable 
membranes that appears during the electric stimulation of tissues by external electrodes 
gives us the clue to define a suitable region of influence, introducing an auxiliary 
concept that in some respects is like the extrapolated distance of nuclear reactor theory. 
From the very beginning it is necessary to consider two different situations: when the 
stimulation is produced by fluid electrodes and when it is produced by a single active 
electrode or by an array of electrodes located near enough to the target element (fibre or 
tissue). 

When we use fluid electrodes, the polarization is usually significant only at the ends of 
the target element or in a region in which the electric properties change rather suddenly. 
The extension of this polarization usually increases with time, and if the applied 
external current field is sustained enough in under threshold conditions, the polarization 
thus produced approaches a steady state. The maximum scope of polarization can be 
estimated using the cable properties of the target element (13). 

A similar situation is encountered when the stimulation is produced by an internal 
microelectrode, as will be briefly considered later (in part I). 

An intermediate case appears in the stimulation of a curved fiber by a fluid electrode, 
introduced in part C. As formula [4d] shows, in this case we can have a non-zero form 
factor in a certain interval of the axis of the fiber and this form factor produces a 
localized polarization of fibre’s membrane similar in that respect to the localized effect 
of a point electrode or array of electrodes which will be considered now. 

When the electrode or array of electrodes is located near enough to the target element, 
the polarization is determined mainly through the spatial variation of the applied 
external electric current field. 

This variation acts on the excitable membranes of the target element through the form 
factor of the applied field: either a volume form factor (as is implied in the right-hand 
members of equations [4a] for a fibre and [5a] for a functional syncytium) or a surface 
form factor (as appears in the right-hand member of the boundary condition [5c] for the 
syncytium). 

The form factor is a scalar function of arc length (for a fibre), of the position in the 
boundary surface (for the functional syncytium) or of the position in a certain three- 
dimensional region (also for the syncytium). Consider the family of points (in the fibre’s 
case), curves (in the boundary of the syncytium) or surfaces (in the interior of the 
syncytium) in which the form factor is constant. 

Next, consider the fibre axis (in fibre’s case) or the lines orthogonal to the curves or 
surfaces (in syncytium’s case), and take the restriction of the form factor to one of these 
lines. We thus obtain the values of the form factor F as a function F..(s) of the arc 
length Scorresponding to the line I’ that we are considering. After that we follow the 
line, leaving the region that is significantly polarized. We will arrive at a final inflection 
point s, of the functionF.(s), after which the form factor will tend monotonically (and 
usually quickly) to zero. Then, we calculate the following estimate of the distance 
(along the orthogonal line I’) in which the absolute value of the form factor suffers a 
significant decrement towards zero: 
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ee IF, (s,) 


F 

ds (s:) 
(This is like the extrapolated distance of nuclear reactor theory). 

Next, we identify the point of I whose arc length is given by s, + As, and we repeat 
the operation for each orthogonal line (in the syncytium case, because in the fiber’s case 
we have only one line: the fibre’s axis). 

Finally, we can define the boundary of the region of influence as the set formed by these 
points. 

Once the region of influence has been constructed, we can apply the procedure 
summarized at the beginning of this section to obtain a coupled system of first order 
nonlinear differential equations in mode amplitudes, as is exemplified in parts E, F and 
G of the present work. 

To identify the threshold manifold in the space of mode amplitudes we must estimate 
the size of a perturbation from the rest state (the origin) which causes instability (in the 
sense that it carries the state of the multivariate system out of the bassin of attraction of 
the origin). 

To be able to make an analytical approach it is necessary to work with a small number 
of dominant mode amplitudes, reducing thus the original multivariate system into a low 
order dynamic. To do this reduction we have at our disposal two procedures. 

The first procedure can be applied when the parameters of our system are close enough 
to a critical value in which the asymptotic stability of the origin disappears. 

Usually, this change of stability is preceded by the critical slowing down of a few mode 
amplitudes (the dominant ones) that slave the others in such a way that (perhaps after a 
short transient) the evolution of the state of the system occurs in a low dimension 
manifold (which exactly at the critical point is the so called centre manifold in 
bifurcation theory). 

When there is only one mode amplitude that shows a critical slowing down, we can use 
Eckhaus method for estimating the size of a perturbation that causes instability. 

In this method the dominant mode amplitude is uncoupled from the others so that its 
threshold behaviour can be studied in isolation. 

Eckhaus method can be expanded and modified to apply it in the framework of the 
second procedure for reducing a multivariate dynamic system to a low order dynamic. 
This second procedure, as the first one, is also framed in a structure of widely separated 
time scales (in our case due to order of magnitude differences in certain kinetic 
parameters that involve some of the membrane’s parameters, the size of the region of 
influence and other geometrical and physical parameters). 

But now, it isn’t necessary to have our dynamic system near a critical point of the rest 
state in the bifurcation diagram, as is needed in the first procedure. In this case the 
equations for the mode amplitudes can be rewritten in a form that allows the application 
of singular perturbation theory in the sense of Tihonov. The fast variables can be 
eliminated in favour of the slow ones, defining thus a slow manifold in which the 
relevant part of the dynamics takes place. (The slow manifold plays the role analogous 
to the centre manifold in bifurcation theory). 

This second procedure is better suited to the study of threshold dynamics by a low order 
modal analysis, because the rest state never loses its stability when the parameters of the 
electrode-target element system vary under normal conditions. 

Some additional remarks can be found in (19) and in the Appendix (where some 
bibliographical references are provided also). 
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5. Aderivation of Lapicque-Hill’s formula for cathodic make excitation of an 
unmyelinated nerve fiber 


The simplest example of the general procedure outlined chapter 4 is afforded by a 
straight and uniform fibre whose membrane can be considered with its recovery 
variables frozen at its rest values and with its activation variables always relaxed to its 
steady-state values. Then, the ionic current density of the membrane is a function of V,, 
alone and can be well approximated by a cubic polynomial. In that case we can work 
with a single nonlinear partial differential equation (Nagumo-Rattay’s equation) (13, 16, 
19) 


ae 


ae 2 ARI, JHA ar Ae +A F(x) 1, (t) [6 a] 
OV Ory. dl, (t) 
ai _ +Ry JVy) =A wy 5 =a eae F(x) [6 b] 


Here x is the coordinate along the axis of the fibre, which we assume to be unbounded 
in both directions. Equation [6 a] corresponds to electric stimulation by an external 
electrode and equation [6 b] corresponds to magnetic stimulation by an external 
stimulating coil. F(x) is the form (or geometric) factor for electrical stimulation (due to 
a single external active electrode) and Fp(x) is the form (or geometric) factor for 
magnetic stimulation (due to a single external stimulating coil or an external array of 
series connected stimulating coils). In(t) is the electric current injected through the 
active electrode and Ig(t) is the current in the stimulating coil. 


If there are several active electrodes (each one with a different current Ip) the forcing 


term in [6 a] is now 4 Tense RNG where the sum extends to the complete set of 
j 


active external electrodes. An analogous expression can be posed for the case of a set of 
external stimulating coils with different known currents through the coils. 


Consider a certain array of stimulating electrodes close to each other and close to the 
fiber. Then the array produces a relatively localized pattern of polarization in fiber’s 
membrane. 


If we introduce the further restriction of considering only one small convex active 
electrode and an infinite length fibre in a homogeneous and isotropic volume conductor, 
then the form factor is given by the following formula that is obtained considering a 
point electrode located above the originx =O, at a distance z, from an unbounded 


straight fiber extended along the axis x : 


1 Oe ag 
5 

. ‘ 9 
4-1-0 (x2 zy y 


F(x) 


When the spatial coordinate tends to +00, the form factor tends as x” to zero. 
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In the case of a compact array of active electrodes and an infinite fiber model, the form 
factors tend to zero as x*. An array of electrodes behaves as a single point electrode far 
enough from the array. The current of this large distance equivalent electrode is the sum 
of the individual currents in the electrodes of the array. 


Before threshold is attained, in an infinite fibre model the perturbation in the membrane 
voltage due to the electrodes behaves as the form factor, tending to zero, when the 
spatial coordinate tends to +00, approximately as x*. From a certain value of |x| 


onwards, the perturbed membrane voltage tends monotonically to its rest value. 
Each monotonic decrement begins at an inflection point of the pattern of membrane’s 
voltage. 


As the perturbation in membrane voltage behaves like the form factor far enough from 
the electrodes, an interval of influence can be defined for a single active electrode, 
following the general procedure shown in chapter 4. When the electrode is near enough 
to the target fiber its size and shape are relevant in the behaviour of the form factor. 


In general for any form factor, an interval |x,, -6,,x;, +5] is obtained, where 


[Fe (x;1.)| [Fe tn) , ; , 
6. = =_——_— , X,, and x,, being the extreme right and left inflection 
 \dF * |dF = oa 
E ; E 2 
ee (XL) dx (Kip) 


points of the given pattern of polarization (Fig.3). 


Fig.3 Sketch of the construction of the interval of influence of an external electrode over a fiber. 
The shape of Fr(x) corresponds to a concave electrode withouth axial symmetry or a symmetric 
concave electrode in a heterogeneous volume conductor. 


Then we assume that the perturbation in the membrane voltage can be zero, at any 
instant, when X = X;p + O,and for x = X,, - 6, 
2. 
If lW, (x)} are the eigenfunctions of the operator — a2 with Dirichlet’s conditions at 
Xx 


the end points of the interval of influence, we can expand membrane voltage as follows: 


Vu (44) = >A, OY, &) [7] 
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If we put R,,J(V,,) = Vy - bVy + CVy, in equation [6], substitute V,, by its 
expansion [7] and project onto the eigenfunctions y,(x) , we obtain an infinite system 
of nonlinear ordinary differential equations for the mode amplitudes A, (t): 


dA hae . 
je er SM )A, +b > Bong Ap Ag - 
dt l p.q=l 


26 Di Ag Re AP Ag By IG) 9 a1 23 0 [8] 


p.g.t =1 


Here Z is the length of the interval of influence of the electrode over the fibre. It is 
nearly equal to 6z (13), where z is the distance between the centre of the electrode I 
and the axis of the fibre. F,,, is the projection of the form factor F,(x) onto the n-eigen- 


function. It is inversely proportional to /* (13). The positive parameters b and c stem 


from the cubic polynomial that gives the ionic current density through fibre’s 
membrane. 


Now, let us consider cathodic make excitation by a small convex electrode in a 
homogeneous and isotropic ohmic volume conductor. Then, as said above, the form 
factor is given by the following formula that is obtained considering a point cathode 
located above the originx=0, at a distance z,from an unbounded straight fiber 


extended alog the axis x : 


1 2x? = 2," 
5 

. ‘i 2 
4-1-0 (x2 zy yp 


Fig.4 shows this geometric factor (or form factor) for a point cathode, the corresponding 
construction of the interval of influence, the Fourier decomposition (the first three 
terms) and its reconstruction (with these terms). 


F(x) 


?7 F(x) r : 
1 
1 
° 1 
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2 x 
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s 4 
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Interval of Influence Fourier’s decomposition Reconstruction of the 
Geometric Factor for a of geometric factor in the | geometric factor in the Interval 
point electrode Interval of Influence of Influence using the first 


three Fourier components 
Fig.4 


Fig.5 shows the results of a digital simulation of the dynamics of the coupled modes 
according to equations [8] truncated after the fifth mode. 
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Fig.5 Digital simulation of the evolution of the first three coupled modes for a fiber stimulated with 
a rectangular current pulse (1ms duration) 


Then the amplitude A, of the first mode happens to be dominant, and as first 
approximation it can be uncoupled from the other mode amplitudes (13, 16, 19, 20). 

With this simplification we obtain a one-dimensional mode amplitude space (the line 
that corresponds toA,) in which the dynamics is given by the following equation (13, 
16): 

Tae 
V’ 


cme ee 


a YA, +b 9, Ay -6%), A, t4y El) [9] 


3, =1,2 and &,,,,=1,5 
If we put /(t) equal to zero, the resulting equation has three equilibrium points: the 
origin (stable), a threshold point A,, (unstable) between the other two, and an excited 
state A,. (stable). In this case the decaying set is given by A, <A,,, the amplifying set 
is given by A, > A,, and the threshold manifold is given by (Fig.6) 


2 
Ao SROs Af La ses) [1+ | [10] 
lu 
2041 2c 41 CA l? 
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Fig.6 A global saddle-node bifurcation of the steady state values of the fundamental mode 
amplitude as a function of the dimensionless quotient (between interval of influence length and 
fiber space constant). The zero-amplitude steady state (rest state) is always present 


Linearizing the equation [9] in the decaying set we obtain: 
2 2 
A 
rai) eae I(t) [11] 


Now, consider a rectangular pulse of amplitude J and duration d that begins with the 
membrane at rest (that is: A,(0) = 0). The condition for a just-threshold pulse is that 


A, (@=A,, [12] 


Integrating [11] for a rectangular current pulse and considering the just-threshold 
condition [12] we finally obtain Lapicque-Hill’s formula: 


a [13] 


an ee 
The rheobase R=|14+2 — ~ [14] 
l Au &, 


Fig.7 and Fig.8 give the rheobase for an unwrapped fiber as a function of electrode-fiber 
distance and as a function of fibres radius. 
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Fig.7 Rheobase as a function of electrode-fiber distance for a fiber without a sheath of connective tissue 
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Fig.8 Rheobase as a function of fibre radius for a fibre without a sheath of connective tissue 


T 
The system’s time constant t= — [15] 
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147 Au 
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Fig.9 and Fig.10 give the time constant for an unwrapped fibre as a function of electrode- 
fibre distance and as a function of fibres radius. 
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Fig.9 Time constant as a function of electrode-fibre distance for a fibre without a sheath of 
connective tissue 


0.00 0.40 0.80 1.20 1.60 2.00 


af 


Fig.10 Time constant as a function of fibre radius for a fibre without a sheath of connective tissue 


Fig.11 and Fig.12 give the limit threshold charge for an unwrapped fiber as a function of 
electrode-fibre distance and as a function of fibres radius. 
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Fig.11 Limit threshold charge as a function of electrode-fibre distance for a fibre without a sheath 
of connective tissue 
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Fig.12 Limit threshold charge as a function of fibre radius for a fibre without a sheath of connective 
tissue 


The fibre model considered in this section can be improved to consider that most fibres 
are wrapped by a purely resistive sheath of connective tissue. The importance of this 
sheath in the distribution of external current through fibre’s membrane grows as the 
electrode gets closer to the fibre. It is possible to derive Lapicque-Hill’s formula also in 
this case. 

As an example of the new results that are obtained, we present the formula for the time 
constant (13): 
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2 
1, =Ty. : [16] 


Here: 4. = fe is a second cable’s space constant (r, the sheath’s resistance per unit 
Tj 


length of the fibre and r; is fibre’s core resistance per unit length). 

Expression [16] can be further improved considering the effects of activation and 
recovery variables. However, such as it is, it can cope with the until now unexplained 
findings of classical electrophysiologists (5) that related the chronaxie with the distance 
from the electrode to the excitable tissue, with electrode radius, with fibre’s radius and 
with the thickness of the resistive sheath that wraps the fibre, amongst others parameters 
(for a convex electrode, /=6- (r, + h), being r, the electrode’s radius and h the 
distance between the electrode’s surface and the target fibre). 

Formula [14] gives the rheobase as a function of the distance between the electrode and 
the fibre (through/) and fiber radius a, (through/,,). From a certain value of / 


onwards, the rheobase R_ increases with /, first linearly and then as 1° (13, 16). 

But when / approaches zero, R increases again, presenting a loop that ends in a point 
with vertical tangent. So, the rheobase, considered as a function of /, presents a 
minimum. 

Now, considering R as a function of the fiber radius a fe it decreases from +0 when 


a, = 0 towards a minimum for a certain value of a pa and then the rheobase increases 


again making a loop that also ends in a point with vertical tangent. 


For a fiber wrapped with a sheath of connective tissue we have another non dimensional 


A : ; . ; 
parameter ——. Here 2, is a second space constant introduced in formula [16] of this 
M 


section. 


A. : ied ; 
If —* is less than a certain critical value, the relation between the rheobase and the 
M 


length of the interval of influence / is qualitatively the same as it is for the unwrapped 
fiber. 


A 
However, when as approaches its critical value from below, the value of / at which 
M 


the minimum occurs and the size of the loop preceding the point with vertical tangent 
both tend to zero. 
A 
If Ss 
Ay 
longitude of the interval of influence | tends to zero 
For the wrapped fiber, the relation between the rheobase R and a, is modified also, in 


is greater tan critical, the rheobase R tends monotonically to zero when the 


AY 


A 
the same manner, depending on the value of : 
M 
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Fig.17 Threshold manifold, decaying set and amplifying set in A1, A3, As space 
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6. Non-linear modal analysis using the Fitzhugh-Nagumo model for the unit 
membrane of the fiber and a derivation of classical two factor models of 
excitation. 


The construction of the impulse response function and the threshold charge for 
electric stimulation, and the impulse response function and the threshold current for the 
magnetic stimulation of a (non-myelinated) fiber will follow a common procedure. This 
allows an easier grasp of the similarities and differences between the two situations. 

Let us begin with the nonlinear cable equation with Rattay’s activating function 
generalized to consider both electric and magnetic stimulation. 

To simplify the derivation, we use the well-known Fitzhugh-Nagumo model for the 
unit membrane. If the (possibly bended) target fiber is represented by a curve of directed 
arc length s measured from a suitable fixed point of the fiber, v(r,s)is membrane voltage 


field and w(t, s) is the recovery variable field, both relative to their rest values, we 
obtain the following set of equations: 


[Sa] 


=(v-w), [5b] 


The time constant of the membrane isz,,, andz,,is the time constant of the recovery 


variable (under physiological conditions at least an order of magnitude greater than r,, ). 
The fiber’s space constant is 2,,. The parameters of the unit membrane model b,c,a,7 are 
positive. 

Rattay’s generalized activating function is given by: 


A, F,,(8)-i, (t) > electric 


7 Ai, Fy (8) < i, (t) > magnetic 


The functions F,(s) and F,,(s) give the spatial distribution of the perturbation produced 
on the excitable membrane by the electric field due to the electrodes or induced by the 
time varying magnetic flux due to the coil, and may be called geometric form factors for 
electric and magnetic stimulation, respectively. Now, the stimulation is always a 
localized phenomenon. Therefore, outside a certain interval of influence along the fiber, 
of length @ , the form factors may be neglected. If we assume that until the fiber reaches 
a threshold state, the fields of membrane voltage and recovery variable take their rest 
values at the end points of this interval of influence, it is possible to make a nonlinear 
modal analysis of the cable equations. The value of ¢ depends of the form factor. If s is 
measured from the mid-point of the interval of influence, we may use the Fourier 
development: 
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+), 
; me [7] 
Fy (s)= Fu | cos Fy2 | sin a 


Substituting [7] in the nonlinear cable equation and eliminating the spatial 
dependence, we obtain a system of first order nonlinear ordinary differential equations 


in the unknown mode amplitudes 4; (1) and B, (t ): Solving this system with suitable 


initial conditions, we obtain the mode amplitudes. The procedure is the same as already 
developed for stimulation by electrodes. The only difference between electric and 
magnetic stimulation in this approach is the forcing term. So, the dynamics of the 
unforced fiber, after the end of the stimulating pulse, is the same if the length of the 
interval of influence is the same. 

Digital simulation of both cathodic make and anodic break electrical stimulation using 
up to seventeen mode amplitudes suggest that the first mode is the most relevant in 
determining threshold behavior of the fiber (40). 

Truncation to the first mode, uncoupled, gives in both cases, a set of nonlinear 
ordinary differential equations in the mode amplitudes corresponding to membrane 
voltage and the recovery variable [19]. 

The study of threshold dynamics done with the afore mentioned equations, with 
parameter’s values within the physiological ranges, shows that a threshold barrier may 


be defined in phase space (A,B; ) such that when the phase point reaches the barrier, 


an action potential emerges. 
The under-threshold behavior can be described with enough accuracy, up to the 
threshold barrier, by the linear system: 


dA dae = 
c Bis (te? SYA, By +y? Fy cl) 
> t 
dB, 
et e(K yi 
ae dt ( Rete i) " 


This approach is the equivalent, for a fiber with a non-uniformly polarized membrane, 
to the classical discussion of Fitz-Hugh for a unit membrane. Fig 18 shows the 
simplified dynamics in the (A,,B,) state space. The results of the digital simulation 


suggest the introduction of a decaying and an amplifying set, separated by a threshold 
curve. 
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Fig 18. Sketch of the dynamics of the system, simplified by the threshold barrier. R is 


the rest state. The decaying set is shown as a shaded area bounded by the double 
arrowed threshold curve. 


The threshold curve is composed by half-straight line taken from the threshold barrier 
and an orbit in the decaying set that is just tangent to the threshold barrier. The 
construction of the threshold barrier can be seen (13). 


The linear dynamic system [8] can be recast in a matrix form: 


Ds S 2 d . >] 

a= AH) Am Fy ict). [9] 

. a2 - 

A, l+a° — — 

Tm Tn 1 
X= = e — 

B . 1 - ' [0 ce 

1 —— se 

WwW Ty 
~1 An 
In the matrix A, the spatial properties appear in the element [108 ° | 


m 


A 
through the non-dimensional number ra The other elements of the matrix A depend 


only on unit membrane properties through the parameters @,/7,7,,.Ty. 


The solution, beginning from the rest state, and until the first arrival to the threshold 
barrier, is: 


2 rd. see 
x(t)= An “Fy {i-(u)-€! 8, )-du [11] 


From the results of digital simulation (summarized in Fig 18) the condition that 
characterizes a just threshold magnetic stimulation is (Fig.5): 
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[12] 


I= Threshold barrier 


Fig.19. Graphical interpretation of Equation [12]. The unit vector Nis normal to the 
threshold barrier. 


There are analytical formulae for the unit vector N and the distance p as functions of 
system’s parameters. 
We put X(t) from [11] into [12] and define: 
(a) The impulse response function: 


ne 
_ 1 

Gy (= [13] 

nN e; 
(b) The threshold current: 
ee: a 
2 >T> 

An i Fy “nN ej 
Then we obtain the just threshold condition for magnetic stimulation in terms of the 
excitation functional: 


I Th,0 — [14] 


t . 
d i(u) 
G,,(t -u)—— du} =I 
Max i M ( ) du Th,0 [15] 
The same procedure applied to electric stimulation gives 
+T tAs> 
Dp ne ey 
QO = G,(t)= ———_ 
din Fea i 6, ce 


So, if the parameters of the unit membrane are the same, the only difference in the 
impulse response function (in the framework of this mathematical model) is due to 
different values of ¢ related with differences in the form factor. 

The activation of peripheral nerves can be studied under the following modeling 
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conditions: the medium can be considered as homogeneous, while the fiber can be 
considered as straight and infinite, but the external electric field varies in a region along 
the fiber. 

In the framework of the present mathematical model, the degree of spatial localization 
of the perturbation of the peripheral nerve membrane due to the external field is given 
by the length of the interval of influence ¢. 

It is possible to show that an approximate formula for the time constant of the 
cathodal strength-duration curve for a nerve fiber is given by: 

T 


m 


(3 = 
8 4] [16] 
l+z = 
0 


The chronaxy is proportional to ts. For a given peripheral nerve, in the electric 
stimulation /is smaller than in the magnetic stimulation case. From [16] it follows that 
chronaxies for electric stimulation should be smaller than chronaxies for magnetic 
stimulation. This explains the experimental findings. 

Equation [16] is obtained from the simplest model for the impulse response function. 
A more realistic model of Gm (given by [13]) is sketched in Fig.20. 


t 


- Time 


Fig.20. Schematic representation of an impulse response function taken from [19]. 


It allows the calculation of the chronaxies both for cathodic and anodic stimulations. 
However, despite the analytical formulae are different, the same theoretical predictions 
about the behavior of chronaxies for electric and magnetic stimulation of peripheral 
nerves are derived from this more complete model. 
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Threshold barrier 


Fig.25 
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7. The region of influence for an electric syncytium stimulated by a point cathode. 


It is possible to make a somewhat similar derivation of strength-duration formulae for 
the case of cathodic stimulation of myocardium by a convex electrode, but the 
mathematical analysis is less simple. (This is so mainly due to the unequal anisotropy of 
cardiac muscle). To simplify both the definition of the region of influence of the 
cathode over the excitable tissue and the nonlinear field equations of bidomain theory, 
let us consider first the following idealized situation. A point electrode is inmersed in an 
isotropic and homogeneous volume conductor, which occupies the whole space. In an 
orthogonal cartesian coordinate system, the electrode is located at x = y = O and 
Z = Z, (Z, positive). The continuum of excitable membranes, which constitute the 
target element of the stimulation, fills the whole half-space given byz < 0. Then the 
boundary of the target tissue is given by the plane z = 0. The continuum of excitable 
membranes is isotropic and homogeneous, so that the field equation [5a] can be 


OV > 
rewritten as follows: pee (oe F: % +J(V,\,W)=CV? Vi, [17] 
f,f,G,G 
C= 1 e 1 e [18] 
f,G, +f, G, 


C is ascalar parameter. Boundary condition [5c] simplifies at its turn to 


EVs 1 
Pe) See eae 19 
Fa ee he a [19] 
Ut 
J (t,x, y,0)= ZoD) [20] 


Am (x? +y~ $257)” 


If p = x’ + y° is the distance to the axis of symmetry of the field, we have a form 


Zo 


[21] 
4m (p? +297)” 


factor F(p)= 


It tends fairly quickly to zero as p_ tends to + 00 (Fig.26). 


F(p) 


Fig.26 Radius of influence of an external cathode on the surface of an excitable syncytium 
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Now we can construct a circle of influence of the electrode over the plane boundary of 
the target tissue, so that in the boundary of the circle we can neglect the polarizing 
effect of the electrode relative to its effect in the interior points of the circle. 

F(p)is a monotonic decreasing function of p with an inflection point when 


Zz 
P=/),= a It diminishes from F(p ;) up to a value than can be neglected, in a 
F(p ;) 2) 
dial dist iven by: ; = 
radial distance given by oF = 620 
dp Pi 


Then the radius R, of the circle of influence can be estimated by (Fig.26): 


F(p ;) 


R,=p,+ ~1,332, [22] 


F 
ip 


Now we must calculate the depth of the polarizing effect in the bulk of the target tissue. 
This depth must be a decreasing function of the distance p of the considered point to 
the z axis. As Fig.27 suggests, we can assume that the region of influence is grossly half 
an ellipsoid of revolution with axesay = b, = Rg, andcy. 


4 


Fig.27 Region of influence due to an external cathode 


Then we only need to estimate c, to define the size of the region of influence. 

As we did it for an excitable fibre, we assume again that the activation variables are 
relaxed to equilibrium with membrane potential and that the recovery variables are 
frozen. Then we can approximate ionic current density by a third-degree polynomial. 
Now we introduce the three-dimensional space constant: 


Ry C 
Ay ={— [23] 
Xm 
Then equation [17] can be rewritten as follows: 
ON 
Tua tN Vw +0Vw =A wo V? Vu [24] 
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Let us suppose that V,, (0, x, y, z) = O (the rest state). Then, because of the cylindrical 
symmetry of this system, V,, will be a function of t, po and z only, and 
Vy (t, Ry, Z) = O for every t and for every Zz. 

In this case: 


OV 0°V 
V? Vy = + (p <8) 4 [25] 
pep Cp OZ 
We can separate the space variables and develop membrane potential as follows: 
Vaults p. 2) = S “A, (tz) v, (p) [26] 
n=1 
The eigenfunctions v,() verify 
i Wd dv 
(p — (p)) =u," V,(p) [27] 
p dp dp 
For every n: 
v, (Ry) =90 [28] 
a 
v,(p)=C, J,(—2) [29] 


0 


J,(¢) is the Bessel function of the first kind and zero order (21) and 
@,<a,<a,<... are the zeroes of the function J,(¢). 
Then: 


o,=— [30] 


The maximum depth of penetration of the polarization of the excitable membranes is 


2 
a, 


Ra. 
and is attained when membrane’s polarization is in a steady-state due to a sustained step 
current. 

Substituting the ansatz A, ,,(z) v,() in equation [24] and linearizing it, we obtain the 
following approximate equation for A, ,,(Z): 


associated with the lowest eigenfunction v,() and its eigenvalue @ Pi = 


da? Ges Sanne 
—, A, ,,(Z) = ss 


A 1 
ay 1 (Z) [31] 


M 


The solution to this last equation, that tends to zero when z > - © is given by 


A,,,(z) = A,,(0)e' where a [32] 


ile aot 
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This parameter can be considered a measure of the maximum depth of penetration of 
membrane’s polarization. So, we obtain: 
Cap ks [33] 


P isa proportionality factor which in a first approximation can be taken equal to one. 


Fig.28 shows , as a function of Re, : 
M M 


O Re 
Ane 


Fig.28 Depth of influence of a point electrode as a function of the radius of influence in the plane 
boundary of an electric syncytium 


The squares of the first three zeroes of the Bessel function are: 


a,” 5.7830 a, ~30.4715 a, ~74.8865 
Only the first is involved in the formula forL, . 
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8. A model of the threshold behavior of electric syncytia. Derivation of strength- 
duration formulas. 


Let now us consider the conditions for the emergence of an action potential in the 
framework of the simplified model of an excitable electric syncytium established in part 
F of this work. 

We begin with an external electrode or array of external electrodes, such that they 
polarize a given bounded region of the continuum of excitable membranes. The 
boundary ¢ B of this region can be partitioned in two regular surfaces: 0 B, which is in 
the bulk of the tissue, and 0 B, which is in the interface between the electric syncytium 
(the bidomain) and the external volume conductor. We assume that V,, = 0 in every 
point of B,. 


ON 
Once the applied current has stopped, we assume that a . = OondB,. 


Here Pn represents the directional derivative in the direction normal to the surface at 
n 


each point of B,. (This homogeneous Neumann’s boundary condition stems from 
equation [5c] written for an isotopic medium when the electric current density vanishes 
in the equivalent volume conductor). 


To obtain an approximation to the amplitude of the localized electric perturbation that 
causes instability and produces an action potential in the tissue, we put 


Vult 1) = 3A, (1) [34] 
n=1 
y ,(1) are the normalized eigenfunctions of the operator — V* with homogeneous 
Dirichlet’s boundary conditions in 0B, and homogeneous Neumann’s boundary 
conditions in 0 B,. 


Let us introduce the inner product (f, g) = Jace) g(r) dV Then, taking into due 
B 


account that (y we WV =) = 6,,, (O,m 1S Kronecker’s delta) and employing the 

; : , ‘ Og Of ; 

identity (22) f Vi gdV = (3 V' fdV + Je — ds - Ds —— dS we substitute 
B B OB 0 n OB 0 n 

[34] in equation [24] and obtain, after several transformations, the following system of 

nonlinear ordinary differential equations for mode amplitudes A , (t) 


d foe) 
ie a fa =U t+Ayo wy, A, +b DT, Ap Ag - 


p.q=l 


ee) 


2 Dili AR We At TO). (BS 


pq@ral 
Here 1,,9=J¥nVpV¥q dV and Tog: =]U¥n Vp Vg Vr AV. 
B B 
The yw, are the eigenvalues corresponding to the eigenfunctions y, 


They verify uw, < uw, < uw; <.. Besides wu,” 
The source terms deserve a special attention: 


—> +0 when n>+0. 
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f, I(t) = ‘| yp, 2% as [36] 
dB, 


If we express the electric current density in the ohmic volume conductor as a function of 
the global electric current produced by the active electrode (or electrodes) we have: 


J(t,r) =I) F(x) [37] 
Here F (1) is a vector field which gives us the spatial distribution of electric current. 
OV, n.d 
As = from [36] and [37] it follows that: 
On EG 
f n. FdS 38 
; "a 7 Ea.’ ! 


Then, : : F(r) , restricted to  B,, can be considered as a “surface form factor” whose 
projection onto yw, gives the direct influence of total applied current over the mode 
amplitude A ,. (There is also an indirect influence of applied current through nonlinear 
mode coupling). 

Let now us consider a convex cathode. In this case we can expect that the amplitude of 
the first mode will be dominant, so that the other mode amplitudes will be slaved by the 
first one. If we further uncoupleA ,, it follows that: 


d 
a Prear =-(1t+rAy we PA,+bI,, A’ -cl,, AP tA wy £10 [39] 


This equation is formaly the same as equation [9] of part (E), now with yw ,” instead of 
2 


—, and with f, instead of F.. 


When I(t) is zero, equation [39] has three equilibrium amplitudes A, = 0, 
bI bia)... a 
A ; = 111 = { 111 oA d 7s y) 7 ul *) 
2cTiy, 2¢cliy Clin eis 


2 
bI bI 1 
eee ia { Wi ) d+ peer ie. 
2Clin 2cTuny Chan 


The first and third one are stable and the middle one is unstable (it is the threshold 
amplitude). The eigenvalue  ,° depends on the shape and size of the region of 
influence of the electrode. When the excitable tissue fills a half-space and the electrode 
is a point or a convex electrode, the region of influence B is half an ellipsoid of 


revolution with a, = b, = R, andc, ~ L, as we showed in part (F). In this case: 
2 > y [24a . 1 [40] 
pe SAE aL Re ee 


(y is anumber between | and 3) (23). 
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For a point electrode Ry, ~ 1,33z,, but this relation can be applied also to a convex 
cathode if z, is the distance between the tissue and the electric centre of the electrode. 
Fig.29 shows the bifurcation diagram as a function of A, “,’. 


Ax, ea 


? 2 
CAPAA), An fa 


Fig.29 Saddle-node bifurcation for the equilibrium values of mode 1 amplitude 


There is a critical value of yz,’ such that once it is attained the threshold amplitude 
disappears and the system cannot produce an action potential. This critical value will be 
always attained with a point electrode if the electrode is close enough to the excitable 
syncytium. When z, tends to +00, threshold amplitude A,, decreases toward a limit 


value. The just threshold condition is max {A,0} = A,, (when the applied current 


te [ 
begins at t = 0). 
Linearizing equation [39] in the decaying set A, < A,,, integrating it from A,(O) = 0 
and using the just-threshold condition for a rectangular pulse of just-threshold amplitude 


I, and duration d, it follows that: Oe [41] 


= (+A ie HL :) fi G, Ai, 


R ——= [42] 
yar Jvin.Fads 
2B, 
i 
fi 43 
Ss tA wy [ ] 


Thus, we have obtained Lapicque-Hill’s formula for the electric stimulation of a 
functional syncytium. Now, if we retain the quadratic term of the evolution equation in 
the decaying set, integrate the resulting differential equation for A, and consider the 
just-threshold condition, we obtain the following threshold relation for a rectangular 
pulse (integrating from A,(0) = 0): 


vote fu $1 [44] 
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Here, by definition: u= Re — |] [45] 
1+Ay “wy f GA 
Now the rheobase is given by R= ( m Hi) fe — [46] 
4blyAw | yw, my. Fads 
OB, 
T M 


The time constant is given by [47] 


t. = ——— 
OA gh 


d 
If we assume that u. s is always much less than 1, from equation [44] we obtain 


Ss 


Lapicque-Weiss formula: x=1l+- [48] 


d 


d 
But this is not a very good approximation for big values of a 


Ss 


I 1 
Now, from [44], after some operations it follows: 7 7 r Fi 7 
2 as 
sen Bi oe | 
d | d ares 
If we put sen k ; a xu. ra and then we find explicitly. we obtain the following 
; I, 1 t- ) 
new formula for the strength-duration relation R* 2! +,/1+4 7 J [49] 


This is a much better approximation to the implicit relation given by formula [44]. 


Now, it is known that if rheobase and limit threshold charge are calculated from 
strength-duration experimental data obtained stimulating cardiac muscle by external 
electrodes, and the resulting predicted curve from Lapicque-Hill’s equation is plotted 
against the experimental results, this curve appears significantly below the experimental 
one at intermediate pulse durations (24). Cardiac muscle membranes have a strongly 
nonlinear current-voltage relation, with inward rectification and very slow recovery (25) 
so that the model that employs a cubic polynomial to express the relation between ionic 
current density and membrane voltage, and even the second-degree approximation used 
to derive formula [44], must be a good approximation in the case of myocardium. Given 
the rheobase and the time constant, the o parameter (introduced in part (A)) is equal to 
1,588 for Lapicque-Hill’s formula [41], equal to 1,618 for the new formula [49] and 
equal to 2 for Lapicque-Weiss’s formula [48]. 

Then, the theoretical curve corresponding to the new formula appears above the 
theoretical curve corresponding to Lapicque-Hill’s formula, but under the theoretical 
curve corresponding to Lapicque-Weiss formula (Fig.30). 
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Fig.30 Comparison between three analytical formulae for the strenght-duration relation: 
Lapicque-Hill, Lapicque-Weiss, and the new formula [49] 


Formula [49] seems to adjust better to the experimental results reported in (24), but the 
whole subject seems to deserve detailed research. Formula [44] can be obtained even for 
an anisotropic and nonhomogeneous domain of excitable membranes (13), but it doesn’t 
take into due account activation and recovery processes. 
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9. Activation and recovery effects: a suggested simple model 


Now we eliminate the restrictions imposed over the activation variables (equilibrium 
with membrane voltage) and recovery variables (frozen). 

Perhaps the easiest way to do this is to introduce one activation variable 7 and one 
recovery variable @ and to assume that their evolution is given by the following linear 
equations: 


s= =H(Vu -9 77) [50a] 


X 


=BWu-7 @) [50b] 


In this approach we suppose that in the rest state 7 = 0, w = 0 and V,, = O. Then in 
equations [4a] for a fiber and [5a] for an electric syncytium, we assume that 
membrane’s ionic current density is given by the following polynomial: 


R,, J(Vy> 1), @) = S10 Vu + 8 77 -aAoO+ 


+214Vm 7 + 8m n° + 21 Vu n° + $03 n° [50c] 


So, the ionic current density depends linearly on the recovery variable w, in the same 
spirit of the well-known Fitzhugh-Nagumo model for a fiber. Moreover, we neglect 
terms of order greater than three in any variable. The parameters a, £, vy, 6 and mw are 


positive and independent of V,,. 

The parameters g,, are related with ionic channel activation, while @ is we related with 
the recovery process. We assume that g,,, g,, and g,, are positive, while g,, and g,, are 
negative. 

Now, we impose to this three-variables model the following restriction: we relax the 
activation variable to equilibrium with membrane voltage, so that 


Then equation [50c] for the ionic current density must reduce to the expression typical 
of Bonhoeffer-van der Pol-Fitzhugh model: 


Ry JVy, @) = Vy - DVy? + cVy? - GO. 


With this restriction, the three variables model proposed in this section can be 
considered as a natural extension of Bonhoeffer-van der Pol-Fitzhugh’s two variables 
model, to take activation effects into account. 

This introduces certain restrictions on the coefficients i (13). We obtain the equations, 


lg 01 | 
ro) 


=1 


IE 10 | ay 
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ul el, 
ro) ro) 

lg. | Eos] _ 
an a 


Please note that here R,, is the electric resistance of the unit area of membrane, in the 
rest state, measured with 7 always relaxed to equilibrium with V,,. The true linear 


membrane resistance R,,, is measured disturbing from the rest state only the membrane 


voltage: aie (0, 0, 0) 
Reine Nag 


R 
From this last definition and from formula [50c] it follows that —“ = Bi 
MO 
Now, we will apply this idealized membrane model to the study of the threshold 
behaviour, both for a single fibre and of an electric syncytium. To do this we must 
substitute, in both cases, equations [4b] and [5b] of section C, by the system of 
equations [50a] and [SOb]. 


We must also substitute R,, J(V,,,W) in equation [4a] (for a fiber) and in equation [5a] 
(for a syncytium), by the simplified formula [50c] introduced in the present section for 
Red as He O)> 

The resulting non-linear field equations can be considered a generalization of the 
Fitzhugh-Nagumo equations for a fibre or for a functional syncytium. 

As was explained in sections C and D and exemplified in sections E, F and G, we pose a 
eigenfunction-eigenvalue problem for a linear self-adjoint and positive elliptic operator 
with pure Dirichlet’s or mixed Dirichlet and Neumann’‘s boundary conditions in the 
region of influence, and we thus obtain a sequence of eigenfunctions y,, each with its 


corresponding eigenvalue yw, . 
Each eigenvalue is inversely proportional to the square of a linear dimension 


characteristic of the given region of influence. 


If the target element is a long fibre stimulated by a relatively close electrode or array of 
2 


electrodes, the elliptic operator is simply - Fre with Dirichlet’s conditions at the end of 
X 


the interval of influence of length 7 . Then, as we saw in sectionE, yw,’ = 


LP’ 


If the target element is a functional syncytium, homogeneous and isotropic - so that a 
single scalar three dimensional space constant is enough - the elliptic operator is - V* 
with mixed boundary conditions (V* is Laplace’s operator) and the eigenvalues are 
given by a generalization of formula [40] for Hie (see section G). When the electric 


syncytium is heterogeneous and/or anisotropic, the elliptic operator is - V. (C.. V) and 


there isn’t a single space constant for the whole tissue. (Some results about this last case 
can be found in (13)). 


In every case the sequence of positive eigenvalues yw, 1s a strictly increasing sequence 


such that, when the characteristic dimension of the region of influence is less than or of 
the same order of magnitude that the fibre’s space constant or the syncytium’s space 
constant (a true space constant for the homogeneous and isotropic case or a suitable 
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defined representative value in the other cases), then each eigenvalue is an order of 
magnitude greater than the precedent one. 
Introducing the ansatz Vi, = DA (Oy, o= > BOY, n= > C,Oy, 


where y, are the eigenfunctions of a Sturm-Liouville problem with homogeneous 
Dirichlet (for a fibre) or mixed Dirichlet and Neumann’s boundary conditions (for a 
functional syncytium), and projecting onto the corresponding eigenfunctions, we obtain 
a nonlinear set of coupled ordinary differential equations for mode amplitudes 
A, @,.B, ©, C,@ @ = 1, 2, 3, «..): 


dA 
Ty ae. ev; (le. - Rie u, )A, 7 Iga C, + lg.,| DP A, C, 
Pq 


Ls Soe os, C, C, 7 Ig. pase A.C, C, 
Pq 


p.q.t 


- [Soa] Do Paar Cp CC, - 2B, +424 FU [51a] 
P.q.t 


dB 
— - Ib 
ae = A (AL 7B.) [510] 
dC 
n = S 1 
Ti u(A,-dC,) [Sic] 


Here the symbols %,, and 3... represent the internal product of yw, with y,y, or 


npqr 
Y.W.,W, vespectively, in an interval (for the fibre’s case) or in a three dimensional 


region (for the syncytium’s case). 


From these equations it is possible to obtain strength-duration formulae which take into 
due account activation and recovery effects, both for nerve or skeletal muscle fiber and 
for myocardium or coupled smooth muscle. Things like the local response, the latency 
in the emergency of the action potential, and the accommodation effects can be thus 
considered from the standpoint of nonlinear modal analysis. Suitable analytical 
formulae can be derived for these phenomena, considering only the mode amplitudes 
A,(t), B,(), C,(t) corresponding to the first eigenfunction, and applying the general 
procedure outlined in section D. Fig.31 shows a sketch of the amplifying set, the 
idealized threshold manifold and the decaying set for a state space constructed from 
mode one amplitudes. 


Fig.31 Threshold manifold, decaying set and amplifying set in Ai, Bi, Ci space 
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Let us consider now why it should be possible to uncouple higher order mode 
amplitudes from low order ones, eventually leaving only mode one amplitudes as a first 
approximation to threshold dynamics. 

We suppose that, before the beginning of the external current pulse, the excitable 
membranes of the target element are at their rest state, so that initially 
A, = B, =C, = 0 forevery n. 

The external current I(t) acts directly on the equations for the mode amplitudesA, , 
belonging to membrane voltage. 

But it acts only indirectly (through the mode amplitudes A.) on the mode amplitudes 
belonging to the activation variables (the C, ) and the recovery variables (the B, ). 
When nincreases the projections F of the form factor of the external electric current 
field tend to zero (except in the case of an internal microelectrode or a synapse). 

The calculations show (20) that it may be possible to have to involve as much as nine or 
ten coefficients F, before they begin to be truly neglectable. 

However, the linear kinetic coefficients in the equation for A, are negative and increase 
very quickly in absolute value when n increases. At the beginning of the applied pulse, 
the source term / ,, F, I(t) begins to shift the mode amplitude A, from its rest value 
(zero). Then the afore mentioned linear term in A, introduces negative feedback that 
tends to stop the evolution of A,. As this braking effect is much stronger on higher 
order mode amplitudes, we can expect that low order mode amplitudes for membrane 
voltage will be significantly more disturbed than higher order ones. Then we can 
uncouple the high order mode amplitudes from the low order ones, and work with the 
later as a first approximation (that is we can work with a low order dynamics in A,, 
B, and C,). 

If we are considering either cathodic make or anodic break excitation of a fibre or 
syncytium with a convex electrode, it is possible to leave only the first order mode 
amplitudes (A,,B,and C, ) neglecting all the others, as a first approach. But if we need 
to study cathodic break or anodic make excitations for a fibre, we need to use at least 
mode amplitudes up to order three to be able to describe what is happening. 

Here, we shall consider that A,,B, and C, are dominant and we will uncouple them for 
all the other mode amplitudes, thus studying the dynamics in a three-dimensional state 
space. We will assume that the higher order mode amplitudes will be slaved by the 
dominant ones. We will not consider renormalization effects due to the influence of the 
adiabatic elimination of higher order mode amplitudes on the above-mentioned 
dominant ones. 

So, we are left with the following relatively simple nonlinear system: 


dC 
ak =p (A,-6 C) [52a] 
dA 
™ a =% (|g.0| +A = Mt, ) A, - lgo| C+ lg.,| Fy, A, C+ [52b] 
Soa! Hi CP ° lg,.| Hin A, C/- Sos] Gin Cr -a B, +A i F I(t) 
dB 
ae =f (A, =f B,) [S2c] 
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Now, let us briefly sketch how we could construct a useful approximation to the 
threshold manifold for equations [52], with the applied external current turned off. 
Under physiological conditions a@/y is big enough, so that the only equilibrium point 


of the dynamic system is the rest state of the membranes: A, = B, = C, = 0 


There is a critical value for ne (that depends on the size of the region of influence in 


the case of a functional syncytium) such that above it there is not an amplifying set. In 
these conditions it isn’t possible to produce an action potential in the target element. 


But if ve is below its critical value, we can construct a good analytical approximation 
to the threshold manifold, and we can establish the conditions for the emergence of an 
action potential. The method for doing this construction will be briefly described. 

: g 1 ; : : 
First notice that the plane C, = 3 A, is a slow manifold that attracts the orbits of the 


dynamic system because C, is a very fast variable in comparison with B, . 


dA 
So, let us freeze B, for a moment. The points of the slow manifold where 4 = 0 


for a given B, are represented by the S curve shown in Fig. 32. 


Z Pat 0, Ci= e A, tine 


Cyt A, pete 
3 


Fig. 32 Slow manifold in Ai, Bi, Ci space 


Between B,, and B,, (with B,, negative and B,, positive) there are three branches of 
that curve, the middle one corresponding to unstable equilibrium values of A,, and the 
others to stable equilibrium values. With B, fixed, the three-dimensional dynamics 
reduces to a two dimensional one in a plane parallel to (A,,C,) . In this plane a 
threshold separatrix divides a decaying set from an amplifying one, as shown in Fig.33 
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Fig.33 Threshold manifold, decaying set and amplifying set in Ai, Ci space 


The two branches of this threshold curve end in the unstable equilibrium 
point (Cy u (B;), Ay u (B,), B, VE 

As B, varies, this threshold point moves on the middle branch of the aforementioned S 
curve, and the corresponding threshold separatrix generates a surface in the three- 
dimensional mode amplitude space. 

If we unfreeze the B, mode amplitude all the equilibrium points 
(C, , (B,), A, (B,)) disappear and only the origin remains as a fixed point of the 
dynamics. Indeed, as was shown in (13), we can approximate this manifold by its 
tangent plane at the point of coordinates C,; = C,,(0), A, = A,,(), and B, = 0. 
This plane can be constructed using closed form analytical formulae and can be 
considered as a kind of excitation barrier that divides the whole state space in two half 
spaces. 

One of these half spaces includes the whole decaying set and part of the amplifying set. 
The other is formed by the amplifying set only. In the first half-space the dynamics can 
be approximated by the linearization of equations [52a], [52b] and [52c] at the origin of 
the state space. 

A cathodic history of external applied current will be just threshold if it leaves the state 
of the system exactly at the reachable subset of the threshold plane. An anodic one will 
be just threshold if it leaves the state of the system at a point such that the future half- 
orbit that begins from this point will be tangent to the threshold plane. 

It’s possible to show (13) that a single threshold condition embraces both situations: 


maximum | G(tt’) I(t’) ar Q.. [53] 


t € [0,+00) 


Here we assume that the external current is applied from t = O on, and that cathodic 
(anodic) currents are negative (positive). 
The Green function G(t) behaves qualitatively as shown in Fig.34. 
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Fig.34 Green’s function (impulse response function) for a model with activation and recovery 
effects 


The first positive maximum corresponds to the completion of the activation process and 
its first negative minimum corresponds to the establishment of the recovery process. 
Q,,< in equation [53] is a certain (negative) cathodic threshold charge. 


Both Q,. and G(t) are obtained as well-defined functions of the parameters of the 


electrode-target element system. 

Equation [53] is a concrete version, derived in the framework of the simple model 
proposed in this section, of the general concept of excitation functional that will be 
mentioned in section I (Final comments and conclusions). This functional summarizes 
the whole threshold behaviour of a target element in the framework of a given 
electrode-tissue system. 
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10. Final comments and conclusions 


The formulae [14] for the rheobase and [15] for the time constant solve the problem 
posed in the introduction, for a single uniform fibre excited by a convex cathode, when 
the activation variables can be considered relaxed to equilibrium with membrane 
potential and the recovery variables can be considered frozen. 


The time constant is approximately proportional to the chronaxie of classical 
electrophysiologists and present-day clinicians. In Fig.35, t, can be seen as a function of 
the distance between the centre of the electrode and the fibre, and in Fig.36, t, is given as 
a function of fibre’s radius, using the improved formula [16]. 


This formula explains the growth of chronaxie with electrode’s radius, its decrement 
with the growth of fibre’s diameter and with the slandering of the sheath of connective 
tissue that wraps the fibre. It also predicts that system’s time constant is always less than 
membrane’s uniform linear time constant 7,, when the activation variables remain 


relaxed to their steady-state values. These predictions seem to be in accord with known 
experimental results. 
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Fig.35 Time constant versus distance between a point cathode and a wrapped fibre, according to 
formula [16] 
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Fig.36 Time constant versus wrapped fibre’s radius, according to formula [16] 
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Even when the improved formulae, which consider the effect of a resistive sheath, seem 
to agree well with the available experimental data, as was mentioned before, a proper 
test of the predictions of the theory outlined here can be made only through thoroughly 
designed new experiments. 


Some preliminary work has been done by Michel Borde and Sebastian Curti, from the 
Department of Neurophysiology, School of Medicine, University of the Republic, 
Uruguay. The results show a good agreement with the theoretical predictions when the 
rheobase, time constant and limit threshold charge are measured in vitro against the 
distance between a point electrode and the median giant nerve fibre of the earthworm. 


Fig.37, Fig.38 and Fig.39 show a sketch of the experimental dispositive and some of 
their results with the cathode close to the fibre. 
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Fig.37 Sketch of Michel Borde’s experimental setup 
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Fig.38 Rheobase as a function of electrode-fibre distance, near the fibre (experimental results of M. 
Borde y S. Curti) 
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Fig.39 Time constant as a function of electrode-fibre distance, near the fibre (experimental results 
of M. Borde y S. Curti). 


Formula [43] for the time constant of an electric syncytium stimulated by an external 
point electrode can be extended to the anisotropic and heterogeneous case. 
Its predictions agree well with known experimental results obtained with cardiac muscle 


(6). 


Formula [42] for the rheobase of an electric syncytium implies that it grows 
approximately as the square of the distance between the centre of the electrode and the 
tissue. But when this distance grows, t, doesn’t tend to r,, for the considered excitable 


tissue but tends to a value which is smaller than z,,. All this also seems to be in accord 
with known and unexplained experimental results. 


Now, there is a point that must be properly stressed: the whole method of nonlinear 
analysis used in this approach to threshold dynamics hangs from the region of influence 
hypothesis introduced in part D. 

This is a significant idealization of what happens in practice because the polarizing 
effect of the external electrodes vanishes gradually, not suddenly as the region of 
influence concept assumes. So, despite the success of its consequences in explaining 
many facts of electric stimulation by external electrodes, a direct validation of this 
assumption, as such, is needed. 


For the case of the interval of influence, this validation was done by Dr. Ruben Budelli 
from the School of Sciences, University of the Republic, by means of a digital 
simulation of the dynamics of a non-myelinated nerve fibre stimulated by en external 
point electrode. 

The obtained results show that the hypothesis of the interval of influence is a good 
approximation when the electrode is not too close to the fibre (20). But even in this last 
case this approximation has not to be necessarily dismissed for the case of before firing 
profiles. 

A similar digital simulation should be done for an electric syncytium stimulated by an 
external point electrode. 
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The main shortcoming of the proposed general procedure is that it is unable to 
cope with the growth of the action potential, being restricted to just-threshold 
dynamics. However, from one side this may be enough for the discussion of 
differential excitation and related subjects of interest from the viewpoint of clinical 
electrophysiology, while from the other side the obtained results could be used as 
guidelines for doing new experiments and new studies of digital simulation of 
threshold dynamics. 


It seems that the idea of an interval of influence can be applied to polarization of fibres 
by internal microelectrodes (13), offering an approach to Rushton’s liminal length (15, 
28) that may be complementary to D. Noble’s one (15) and allowing a theoretical 
analysis of the relation between chronaxie determined by intracellular and by external 
electrodes. 

Using Denis Noble’s model of the steady threshold profile of polarization of a fibre by 
an internal microelectrode, and applying to this profile the same idea outlined in section 
D to define an interval of influence from the form factor, we obtain for its length the 
expression / = 2(A,, + X,) | where 2x, is the length of the segment of excitable 
fibre which is above membrane’s uniform threshold voltage. 

This parameter x, can be related with excitable membrane’s parameters, for example as 
has been done by Noble (15). Then, formula [15] can be applied substituting Zby 
2(A y, + X,). When a point electrode is far enough from the fibre, Z will be greater 
than 2(/ ,, + X,,) so that in this case t, for extracellular stimulation should be greater 
than t, for stimulation by an internal microelectrode. But if the point electrode gets 
closer to the fibre, t, will first diminish, be equal to t, for internal stimulation and, for 
the electrode even closer, it will be even smaller than the time constant obtained with 
the internal microelectrode. This theoretical prediction seems to explain some of the 
known and unexplained experimental results reviewed by Ranck (27). 


The method outlined here can be applied, in principle, to the excitation of nerve, 
skeletal muscle, myocardium and smooth muscle, not only when the active electrode is 
convex but also when it is concave. 

It is possible to determine theoretically the form factor produced by a concave electrode 
in each fiber (13), or its equivalent in a functional syncytium (29). Then, the observed 
difference in behaviour between strength-duration curves produced by convex and 
concave electrodes used to stimulate the same myocardium (30) can be at least 
qualitatively explained. 


In fact, we have been considering tissues stimulated by an active external point 
electrode. But nonlinear modal analysis of threshold dynamics can be applied in other 
cases too. Indeed, when the target excitable tissue is near enough to electrode’s surface, 
electrode’s shape may have an influence on the global parameters of the measured 
strength-duration curves. 

In the framework of nonlinear modal analysis, this influence appears only through the 
projections of the form factor onto the mode’s eigenspaces. 


Now, considered as a distributed electric current source, it seems that any finite 
electrode can be well represented by the first three or four terms of a multipolar 
expansion (13, 29) (1.e., monopole, quadrupole and so on, because the dipole term 
vanishes when the origin is chosen to be at the electric center of the distribution of 
electric current sources). 
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Because of this, each projection of the form factor can be expressed as a convenient 
addition of a monopolar term, a quadrupolar term, and so on. 

If the electrode’s head is concave, instead of convex, and if the target fibre is located 
near the electrode’s surface, the influence of quadrupolar and higher order multipolar 
terms on the form factor and its projections cannot be neglected because they produce 
an important qualitative change in the form factor. 

Indeed, when the target fibre is orthogonal to the electrode’s axis of symmetry, the 
central depolarizing interval typical of a point cathode (or, more generally, of a convex 
cathode) is substituted now by a central hyperpolarizing interval between two 
symmetric depolarizing ones. 

Then, even in the case of cathodic make excitation, thresholds dynamics must be 
studied using higher order mode amplitudes. 

Of course, if the fibre is withdrawn far enough from the concave electrode, the effect of 
the quadrupolar and higher order multipolar terms vanish, and mode one amplitude 
returns to be dominant. 

For a convex electrode the form factor is always like the form factor of a point 
electrode, because in this case quadrupolar and higher multipolar terms don’t produce 
any qualitative variations. 


An excitation functional has been introduced recently (12, 31) to state a just-threshold 
condition that may be used in the study of an optimal pulse shape for pacing excitable 
tissues with external electrodes. 

With the proposed method for the nonlinear modal analysis of threshold dynamics, it is 
possible to derive analytical formulae for the Green’s function (impulse response 
function) of the afore mentioned functional (13, 16, 20). 


Equation [4a] must be suitably modified to be applied to a fibre wrapped with a highly 
resistive sheath of connective tissue (13). Even with membrane activation variables 
relaxed and membrane recovery variables frozen, the threshold behaviour of the 
wrapped fibre shows interesting differences with the situation considered in part E of 
this paper (32). 


Returning to the two basic questions of clinical stimulation stated in section A of this 
work, note that we could plot in the same parameter plane the predicted strength- 
duration curves for different target fibres located in the same electrode tissue system. In 
principle we could determine which elements will be excited by a given applied pulse of 
external current of given amplitude and duration. If two of these curves cross each 
other, we can obtain a differential stimulation of one of the corresponding two elements 
using the fact that one of these elements is, as a part of the given electrode-target 
element system, more excitable than the other at short durations, and the opposite 
happens at large durations of the stimulus. Besides this, the present theoretical 
framework allows the discussion of the possibility of employing suitable histories of 
applied current or modifying the array or type of electrodes to produce the differential 
stimulation of certain target elements in each tissue. 


Finally, the concept of family of threshold states of a given target element (fibre or 
functional syncytium) can be used outside the framework of electric stimulation by 
external or internal electrodes. In fact, it could be applied to the study of synaptic 
excitations, pacemaker activity, and ectopic depolarizations in excitable cells or tissues. 
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Synaptic junctions can be considered much like internal microelectrodes, so that for 
them the length of the interval of influence, when membrane’s activation variables are 
relaxed to equilibrium with membrane voltage and recovery variables are frozen, should 
be 2(4 yy + X,) 


This idea could be applied to study the emergence of an action potential in a single 
skeletal muscle fibre, stimulated by a motor neuron through an end plate neuromuscular 
junction. 


The manifold electric currents produced by the transmitter-gated channels of the local 
muscle fibre membrane could be lumped in a single global current I(t), considered as 
given and injected in a single point inside the fibre. The corresponding form factor 
would be proportional to a Dirac’s delta, cantered at the point of current’s injection. 
Then, a nonlinear modal analysis of threshold dynamics could be done, following the 
guidelines stated in this work. 


The integration of excitatory and inhibitory synaptic signals by a central neuron with the 
emergence of the action potential in the axonic cone is outside of the scope of the 
interval of influence concept. However, an interval of influence may be defined for the 
axonic cone. In any case, the general concept of family of threshold states could be 
applied to this case also. 


For three-dimensional electrically coupled cell networks, each threshold state involves a 
certain space region of excitable tissue that must be depolarized above the uniform 
membrane threshold of voltage if an action potential is to be produced. 


The possible sizes and shapes of these critical regions can be calculated by an extension 
of the analysis given in section G of this paper. 


This would allow us to discuss ectopic excitations and physiological pacemaker from a 
different standpoint, using the concept of critical region to define suitable security 
factors and to understand (in part, at least) the interplay of manifold physiological 
variables that have an influence on threshold. 
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11. Appendix 


(a) The extrapolated distance mentioned in part D of the present work, as well as other 
concepts of the physics of nuclear reactors which have certain interest for a comparison 
with the threshold behaviour of excitable tissues, can be found in any good book about 
nuclear engineering. See, for example, the book by Duderstadt and Hamilton (33). 


(b) The original method of Eckhaus for estimating the size of a perturbation which 
causes instability can be found in (34) and (35). The two general procedures that ensure 
a reduction of a multivariate system into a low-order dynamics - the use of the center 
manifold theory to identify the critical slowing down of a few modes and the use of the 
slow manifold theory to identify a few dominant modes far from criticality - can be 
found in the book of Nicolis (36). Tikhonov’s approach to singular perturbation theory 
giving a foundation to the slow manifold method can be found in (37) and (38). 


(c) It may be of certain interest to compare threshold dynamics with nuclear reactors 
dynamics, because of noticeable similitudes and differences between them. 


It is well known that, for a reactor core of a given composition, a given shape and a 
given intrinsic multiplication factor greater than 1, the core’s dimensions must be 
greater than certain critical ones. If they are not, the nuclear chain reactions cannot be 
sustained (the leakage of neutrons by diffusion will be dominant) and the neutron 
density flux will gradually vanish, no matter how big the initial value of the neutron 
density is. 


In the case of an excitable tissue, we have seen that if the dimensions of the region of 
influence are under certain critical values, membrane voltage will tend to rest without 
producing an action potential, no matter how big is the initial depolarization (for 
cathodic make excitation) or hyperpolarization (for anodic break excitation. 


But when both, the region of influence of the applied electric current field and the 
extrapolated nuclear reactor core, are above their critical size, a significant difference 
appears. 


The zero-neutron density now becomes an unstable equilibrium state. Because of this, if 
we impose zero neutron density at the boundary of the extrapolated core, the neutron 
density flux will always tend to a certain stable steady profile, no matter how small the 
initial distribution of neutrons may be: this is the so called “hair trigger effect”, well 
studied in the framework of Fisher’s reaction-diffusion equation (39). 


When the region of influence is supercritical, the rest state of the excitable membrane 
remains stable (it doesn’t suffer a change of stability as happens with the zero-neutron 
flux in the reactor’s case). 


Now we have a threshold value - the well-known uniform threshold of membrane 
voltage - which is the unstable equilibrium point of the homogeneous kinetic equations 
for the membranes (that is, the equations that remain when spatial differences are 
neglected). 
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Because of this, if a cathodic excitation is going to be successful, the initial pattern of 
membrane voltage has to be well above this uniform threshold in a whole region whose 
dimensions must be greater than certain critical ones. So, we find for excitable tissues a 
second kind of critical region that generalizes the concept of liminal length due to 
Rushton. 


This opens an interesting road for research, in the framework of the concept of family of 
threshold states, introduced in this memory. 


(d) The hypothesis about the dominance of the amplitude of mode one has been verified 
by H. Korenko, of the Faculty of Sciences, University of the Republic, by digital 
simulation of cathodic and anodic stimulation of nerve fibres by an external point 
electrode (40). For a cathodic stimulation of a wrapped fibre it was verified by G. 
Artucio and J. Baraibar, of the Faculty of Engineering, Catholic University of Uruguay 
(32). 


(e) The results obtained in (32) for a wrapped fibre differ from the results of section E 
of this paper only when the electrode is very close to the fibre. 


(f) 


Fig.40 shows some results of the digital simulations done by R. Budelli, for a point 
electrode located at 20 yam and 100 um, respectively, from a given fibre. 


He used the program Neuron, developed by Dr. Hines, to simulate the behaviour of a 
nerve fibre stimulated by a point electrode located at several distances from the fibre. 
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Fig.40 (a) Action potential when the distance between electrode and fibre is 20 Lan 


(b) Action potential when the distance between electrode and fibre is 100 Lan 
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When the point electrode is located closer than 20 am, the interval of influence doesn’t 


corresponds very well to that constructed from the form factor of the activating 
function. 


As it length doesn’t tend to zero with the distance between the electrode and the fibre, 
we shouldn‘t expect to find in true experiments the singular points whit vertical tangents 
shown in Fig.7 and Fig.8. 


Similarly, we can’ t expect that the method for the construction of a region of influence 
of a point electrode over an electric syncytium (like myocardium) proposed in section F, 
works properly when the electrode is located too close to the surface of the excitable 
tissue. This is a subject that deserves careful research. When the electrode is close 
enough to the target tissue, the size and shape of both, the electrode, and the excitable 


61 


Research Report February 1997 


cells, must be considered. As a consequence it seems that the suggested research should 
be done outside the framework of classical one-dimensional cable theory and its naive 
extension in bidomain theory. 
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